clc;
clear;

tic;
startTime = cputime;
pde = process_pde_2d(1);
basis_type = "P1";
disp(["basis_type: ", basis_type]);
vector_h = [1/8 1/16 1/32 1/64];

error_Linf = zeros(size(vector_h,2),1);
error_L2 = zeros(size(vector_h,2),1);
error_H1 = zeros(size(vector_h,2),1);

for i = 1:size(vector_h,2)
    h = vector_h(i);
    mesh = process_mesh_2d(pde, h, basis_type);

    A1 = assemble_stiffness_matrix_2d(mesh, "dx", "dx", pde.lambda);
    A2 = assemble_stiffness_matrix_2d(mesh, "dx", "dx", pde.mu);
    A3 = assemble_stiffness_matrix_2d(mesh, "dy", "dy", pde.mu);
    A4 = assemble_stiffness_matrix_2d(mesh, "dy", "dx", pde.lambda);
    A5 = assemble_stiffness_matrix_2d(mesh, "dx", "dy", pde.mu);
    A6 = assemble_stiffness_matrix_2d(mesh, "dx", "dy", pde.lambda);
    A7 = assemble_stiffness_matrix_2d(mesh, "dy", "dx", pde.mu);
    A8 = assemble_stiffness_matrix_2d(mesh, "dy", "dy", pde.lambda);
    A = [A1+2*A2+A3, A4+A5;
        A6+A7, A8+2*A3+A2];
    [b1, b2] = assemble_load_vector_2d(mesh, pde);
    b = [b1; b2];
    [A, b] = process_boundary_2d(A, b, mesh, pde);
    U = A\b;

    error_Linf(i,1) = process_error_2d(mesh, pde, U, "Linf");
    error_L2(i,1) = process_error_2d(mesh, pde, U, "L2");
    error_H1(i,1) = process_error_2d(mesh, pde, U, "H1");

    fprintf("完成\t%d/%d\n", i, size(vector_h,2));
end

process_convergence_h_2d(error_Linf, vector_h);
process_convergence_h_2d(error_L2, vector_h);
process_convergence_h_2d(error_H1, vector_h);

elapsedTOCTime = toc;
elapsedCPUTime = cputime - startTime;
disp(["TOC time(s)", num2str(elapsedTOCTime)]);
disp(["CPU time(s)", num2str(elapsedCPUTime)]);

%% assemble_stiffness_matrix_2d
function matrix_A = assemble_stiffness_matrix_2d(mesh, D1_type, D2_type, cor)
matrix_A = sparse(mesh.Nb,mesh.Nb);
for n = 1:mesh.Ne
    E = mesh.node(mesh.elem(n,:),:);
    V = E(1:3,:);
    [gauss_points_loc, gauss_weights_loc] = gauss_lengdre_loc(V, "2d");
    phi1 = basis_function_loc(V, gauss_points_loc(:,1), gauss_points_loc(:,2), mesh.basis_type, D1_type);
    phi2 = basis_function_loc(V, gauss_points_loc(:,1), gauss_points_loc(:,2), mesh.basis_type, D2_type);
    for i = 1:mesh.Nlb
        for j = 1:mesh.Nlb
            f = cor * phi1{j} .* phi2{i};
            r = gauss_weights_loc * f;
            matrix_A(mesh.elem(n,j), mesh.elem(n,i)) = matrix_A(mesh.elem(n,j), mesh.elem(n,i)) + r;
        end
    end
end
end

%% assemble_load_vector_2d
function [vector_b1, vector_b2] = assemble_load_vector_2d(mesh, pde)
vector_b1 = zeros(mesh.Nb,1);
vector_b2 = vector_b1;
for n = 1:mesh.Ne
    E = mesh.node(mesh.elem(n,:),:);
    V = E(1:3,:);
    [gauss_points_loc, gauss_weights_loc] = gauss_lengdre_loc(V, "2d");
    phi = basis_function_loc(V, gauss_points_loc(:,1), gauss_points_loc(:,2), mesh.basis_type);
    ff1 = pde.f1(gauss_points_loc(:,1), gauss_points_loc(:,2));
    ff2 = pde.f2(gauss_points_loc(:,1), gauss_points_loc(:,2));
    for j = 1:mesh.Nlb
        f1 = ff1 .* phi{j};
        f2 = ff2 .* phi{j};
        r1 = gauss_weights_loc * f1;
        r2 = gauss_weights_loc * f2;
        vector_b1(mesh.elem(n,j),1) = vector_b1(mesh.elem(n,j),1) + r1;
        vector_b2(mesh.elem(n,j),1) = vector_b2(mesh.elem(n,j),1) + r2;
    end
end
end